Efficient simulation of infinite tree tensor network states on the Bethe lattice 
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We show that the simple update approach proposed by Jiang et al [H.C. Jiang, Z.Y. Weng, and T. Xiang, 
Phys. Rev. Lett. 101, 090603 (2008)] is an efficient and accurate method for determining the infinite tree 
tensor network states on the Bethe lattice. Ground state properties of the quantum transverse Ising model and 
the Heisenberg XXZ model on the Bethe lattice are studied. The transverse Ising model is found to undergo a 
second-order quantum phase transition with a diverging magnetic susceptibility but a finite correlation length 
which is upper-bounded by 1/ \n{q - 1) even at the transition point {q is the coordinate number of the Bethe 
lattice). An intuitive explanation on this peculiar "critical" phenomenon is given. The XXZ model on the Bethe 
lattice undergoes a first-order quantum phase transition at the isotropic point. Furthermore, the simple update 
scheme is found to be related with the Bethe approximation. Finally, by applying the simple update to various 
tree tensor clusters, we can obtain rather nice and scalable approximations for two-dimensional lattices. 
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I. INTRODUCTION 

The investigation of quantum lattice models is one of the 
central themes in modern condensed matter physics. A cru- 
cial step is to develop novel numerical methods to efficiently 
simulate the interesting and complex phenomena of quan- 
tum many-body systems. In particular, the tensor network 
states and the related renormalization group methods, includ- 
ing the tree tensor network state (TTN),'~ 3 the multi-scale en- 
tanglement renormalization ansatz, 4 the projected entangled 
pair state, 5 the tensor renormalization group (TRG), 6-8 and 
the second renormalization group (SRG), 910 are now under 
rapid development. These methods provide promising numer- 
ical tools for studying strongly correlated systems, especially 
for the frustrated magnetic systems and fermion models, and 
can be regarded as an extension of the fruitful density ma- 
trix renormalization group (DMRG) 1 1 in two or higher dimen- 
sions. 

In the study of the tensor network methods, one needs 
to first determine the tensor network wavefunction for the 
ground state. In Refs. [7,10], a simple update scheme is pro- 
posed to determine the ground state tensor network wavefunc- 
tion in two dimensions. This scheme is efficient and robust. 
It proceeds in three steps: (1) apply the imaginary time pro- 
jection operators simultaneously on bonds of the same type, 
for example the x directional bonds in Fig. la, and enlarge 
the bond dimension; (2) construct a local evolving block ma- 
trix and simulate the environment contribution by the diag- 
onal matrices on the external bonds [A y and y z in Fig. lb]; 
(3) decompose the evolving block matrix by singular value 
decomposition (SVD) and decimate the vector space of the 
enlarged geometric bond according to the singular values in 
the updated diagonal matrix 6\. This technique has been com- 
bined with the TRG/SRG to evaluate the ground state prop- 
erties of two-dimensional (2D) Heisenberg models. 8 ' 10 ' 12 ~ 14 It 
is an accurate numerical method for evaluating local physi- 
cal quantities, but it is less accurate in evaluating the long- 
range correlation functions. 10 This is the major drawback of 



this simple update scheme. It results from a mean-field ap- 
proximation for the environment tensor. A way to go beyond 
this approximation is to enlarge the size of the cluster that is 
used for evaluating the environment tensor. This, as shown by 
Wang and Verstraete 15 , can indeed improve the accuracy for 
the long range correlation function. 

In this work, we apply the simple update scheme to infinite 
TTN (iTTN) states on the Bethe lattice. We will show that 
this is a quasi-canonical approach for treating an iTTN. Here 
by the word "quasi-canonical" we mean that with increasing 
the number of iteration steps and decreasing the Trotter error, 
the tree tensor network state obtained by the simple update 
scheme would become asymptotically canonical [i.e. the ten- 
sors satisfy certain canonical orthonormality conditions, see 
Eqs. (3) below]. Thus the simple update scheme provides an 
accurate and efficient approach for evaluating the ground state 
wavefunction on the Bethe lattice. 

The Bethe lattice, as shown in Fig. la, has a self-similar 
structure with an infinite Hausdorff dimension. The size of 
the lattice is infinite, hence the boundary effects do not need 
to be explicitly considered. The Bethe lattice was first used in 
the study of classical statistical mechanics. 16-18 It has attracted 
broader interest since a number of chemical compounds with 
the Bethe lattice structures, such as the dendrimers, 19 have 
been synthesized in the laboratory. 20 

A finite Bethe lattice is called a Cayley tree. Soon af- 
ter White's invention of DMRG, 11 the DMRG algorithm for 
the quantum lattice models defined on the Cayley tree was 
proposed. 21,22 Based on the DMRG calculation of local phys- 
ical quantities in the central part of the Cayley tree, Otsuka 21 
claimed that the anisotropic S=l/2 Heisenberg model (i.e. the 
XXZ model) on the Bethe lattice should exhibit a first-order 
phase transition at the isotropic point. Later Friedman 22 pro- 
posed an improved DMRG scheme and evaluated the spin- 
spin correlations in the ground state. Based on the DMRG 
result, he suggested that long-range magnetic order might ex- 
ist at the isotropic Heisenberg point. Recently, Kumar et al. 
calculated the magnetization with a further improved DMRG 



algorithm, and showed that such long-range magnetic order 
does exist at that point. 

The above DMRG calculations were done on the Cayley 
tree lattice, not on the true infinite Bethe lattice. Furthermore, 
it should be pointed out that the boundary effect is very strong 
on a finite Cayley tree since more than one-half of the total 
sites reside on the lattice edge. This may strongly affect the 
properties of the system. In some cases, the results obtained 
on a Cayley tree lattice can be completely different from those 
for the corresponding Bethe lattice. For example, the classical 
Ising model shows a phase transition on the Bethe lattice, but 
not on the Cayley tree lattice. 24 

To unambiguously resolve the above problems, it is neces- 
sary to calculate the spin models directly on the Bethe lattice. 
The recent development of the TTN algorithms 3 has indeed 
made this feasible. 25,26 In particular, Nagaj et al in Ref. 25 
extended the infinite time evolving-block decimation 27 tech- 
nique to the Bethe lattice and determined the ground state 
wavefunction by imaginary time evolution. For the transverse 
Ising model on the Bethe lattice, it was found that a second- 
order quantum phase transition exists at a critical transverse 
field. An interesting result revealed in this calculation is that 
even at the second-order critical point, the correlation length 
remains finite. For the Bethe lattice with coordination q — 3, 
the correlation length is shown to be less than l/ln2. How- 
ever, in the calculation Nagaj et al used a three-site projection 
operator to simultaneously evolve the two equivalent incom- 
ing legs of the tensors, the computational cost is thus very 
high. The computational time scales as 0(D 8 ) with D the ten- 
sor dimension, which limits the value of D that can be handled 
to D < 8. 

Recently, Nagy 26 proposed a different algorithm to reduce 
the computational cost by making use of the Cj, rotational 
symmetry of q — 3 Bethe lattice. This algorithm reduces 
the computational cost to 0(D 4 ) hence greatly improves the 
efficiency. It can be used for studying the spin- 1/2 quantum 
lattice models. However, the application of this method is re- 
stricted to the translation invariant spin- 1/2 system. 

As will be shown below, the simple update scheme is very 
efficient. Its computational costs scale as O(Z) 4 ), similar as 
for the algorithm proposed by Nagy 26 . But it is much more 
flexible. It can be applied to treat arbitrary TTN states, with or 
without translation invariance. Here we studied two spin mod- 
els defined on the Bethe lattice. One is the transverse Ising 
model and the other is the antiferromagnetic XXZ Heisenberg 
model. The quantum phase transitions and the ground state 
phase diagrams of these models are studied. 

The rest of the paper is arranged as follows. An introduc- 
tion to the simple update scheme and its relationship with the 
Bethe approximation is presented in Sec. II. The study of 
the quantum phase transitions of transverse Ising and XXZ 
Heisenberg models are presented in Sees. Ill and IV, respec- 
tively. In Sec. V, the present scheme is generalized to larger 
tree tensor clusters, in order to provide more accurate approx- 
imations for 2D lattices. Finally, Sec VI is devoted to a sum- 
mary. 
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FIG. 1: (Color online) (a) The q = 3 Bethe lattice. Every site has 
3 nearest neighbors, and the three bonds are labeled according to 
their directions as x, y, and z, respectively, (b) The two-site cluster 
used in the single-bond projection of the simple update scheme. The 
diagonal matrices A y and y z on the dangling bonds should be included 
in the projection to mimic the entanglement renormalization of the 
environment to this two-site system, (c) A minimum cluster that is 
used in the Bethe approximation. It consists of one A tensor and 
three nearest-neighbor B tensors (or vice versa). 



II. THE CANONICAL FORM AND THE SIMPLE UPDATE 
SCHEME 

The iTTN state on the Bethe lattice comprises 4-indexed 
tensors A™ and B™ y defined on the vertices, and the di- 
agonal matrices 9, A, y defined on the bonds along the x, y, 
and z directions as shown in Fig. la, respectively. The bond 
indices represent the quantum numbers of the virtual basis 
states. The physical index m runs over the d basis states of 
the local Hilbert space at each lattice site. The diagonal matri- 
ces store the entanglement information, and play an important 
role in the simple update scheme. 

In order to determine the ground state wavefunction, the 
imaginary time evolving operators U(r) = exp (-r/jy) are ap- 
plied to the iTTN iteratively. At each step, the dimension of 
the evolved bond is increased by a factor of d 2 . Thus the 
tensor dimensions will proliferate exponentially with the in- 
creasing number of projection steps. In order to sustain the 
projections until the iTTN converges to the true ground state 
wavefunction, one needs to truncate the bond dimension after 
each projection step. This needs a proper consideration of the 
renormalization effect of the environment tensor. 

An accurate and full determination of the environment ten- 
sor is computationally costly. This limits generally the tensor 
dimension D that can be handled to a rather small value, say 
D < 6. The simple update scheme 7 ' 9 ' 10 , on the other hand, 
takes the product of the dangling bond matrices as a mean 
field approximation to the environment tensor. It converts 
the complicated global optimization problem into a local one, 
and hence greatly simplifies the calculation. On the regular 
2D lattice, the bond matrix is an approximate measure of the 



3 



entanglement between the system and environment tensors. 
However, as will be shown later, the square of the diagonal 
bond matrix on the Bethe lattice is the eigenvalue of the re- 
duced density matrix if the iTTN is canonicalized, i.e., if the 
tensors in the network are always kept in canonical form by 
some transformations. Thus the simple update scheme is an 
accurate treatment for the renormalization of the iTTN on the 
Bethe lattice. 

The simple update scheme is also closely related to the fa- 
mous Bethe approximation. 16-18 To understand this, we show 
in Fig. lc a 4-site cluster, which contains one A tensor and 
three B tensors. In the simple update calculation, the two lo- 
cal tensors, A and B, and the three inner bond matrices (6 X , 
A y , y z ) should be evaluated and updated iteratively. After each 
single projection step on the inner bonds, to keep the scheme 
self-consistent, one should also update all the dangling bonds 
of the cluster, by replacing the bond matrices with the corre- 
sponding ones on the inner bonds. 

This cluster structure and the self-consistent scheme is in 
fact the Bethe approximation that was first proposed by Bethe 
in 1930's, in the context of statistical mechanics. 16 The key 
idea is to treat the correlations between the central spin and its 
nearest neighbors in the cluster exactly, and to use an effective 
mean field to approximate the interactions between the cluster 
and the rest lattice spins. By solving this simple cluster prob- 
lem, and assuming that all the spins in the cluster have exactly 
the same local magnetization, one can determine the sponta- 
neous magnetization self-consistently. For the quantum cases, 
the six diagonal matrices on the dangling bonds of the cluster 
are taken as the mean fields acting on the inner block. The 
self-consistent condition requires that the matrices 6, A, and y 
on dangling bonds are equal to the corresponding matrices on 
the inner bonds between A and B tensors (see Fig. lc). 

A tensor network state contains redundant gauge degrees of 
freedom on each bond. It is invariant if one inserts a product 
of two reciprocal matrices on a bond and absorbs separately 
each of them to a local tensor at the two ends of the bond. 
This gauge invariance of a tensor network state can be used 
to simplify the calculation of local tensors, especially for the 
iTTN states on the Bethe lattice, where a special gauge, called 
canonical form, can be introduced. 

To be specific, the local tensors of canonical iTTN states 
satisfy the following orthonormality conditions 



m x,y 



(i) 

(2) 
(3) 



where T represents the A or B tensor. If we cut an arbitrary 
bond to divide the Bethe lattice into two parts, denoted as a 
system and an environment subblock, one can then define the 
reduced density matrix of the system block by integrating out 
all the degrees of freedom of the environment block. For the 
tensors that satisfy Eqs. (1-3), the square of the diagonal bond 
matrices are the eigenvalues, and the renormalized bond bases 
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FIG. 2: (Color online) One iteration step in the simple update 
scheme: (a) A and B tensors are connected by the bond x, which will 
be involved in the following projection steps. There are diagonal ma- 
trices on the dangling bonds (y and z bonds) of A and B tensors, (b) 
Absorb the four dangling matrices into A and B, and define the block 
matrix M^py Then take the QR decomposition for M a (b), obtaining 
Q a (b) and R a ^ matrices, (c) Project U(f) onto the bond by contrac- 
tions, and obtain the block matrix G [see Eq. (7)]. (d) Take singular 
value decomposition of G to find the unitary matrices U and V T , and 
the new diagonal matrix ff . (e) Truncate the x-bond dimension to D 
according to the diagonal values of 9'. Merge U (V T ), y' 1 , and Az [ 
together into Q a (b), finally we arrive at the updated A'(B') tensors. 



are the eigenvectors of the corresponding reduced density ma- 
trix, which are orthonormal to each other. Thus, in terms of 
the Schmidt decomposition, the square of the diagonal ma- 
trix elements represent the probability amplitudes of the cor- 
responding eigenvectors appearing in the wavefunction. 

The existence of this simple canonical form of the iTTN, 
i.e. Eqs. (1-3), is very useful in the calculations. First, the di- 
agonal bond matrix describes the entanglement spectrum be- 
tween the system and environment subblocks. Thus to select 
the virtual bond basis states according to the values of these 
diagonal matrix elements provides an optimal scheme to trun- 
cate the bond dimension. Second, the contribution of the en- 
vironment tensors can be faithfully represented by the 4 diag- 
onal matrices on the dangling bonds surrounding the central 
bond under projection (see Fig. 2a). It means that the imagi- 
nary time evolution on each bond can be done rigorously and 
locally. Furthermore, we can also evaluate the expectation 
value of a local operator simply by contracting a small cluster 
comprising those tensors and bond matrices on which the op- 
erator acts. This significantly reduces the computational cost. 

Bearing in mind the benefits of the canonical iTTN states, 
one can perform explicitly the canonical transformations dur- 
ing the projection processes. However, to further save com- 
putational costs, in practical calculations we choose to carry 
out the imaginary time evolution just using the simple update 
scheme, and gradually reduce the Trotter step r, which would 
bring the iTTN states into its canonical form step by step. This 
scheme works because the diagonal bond matrix provides an 
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approximate measure for the entanglement between the two 
sides of the bond and can be used to substitute approximately 
the environment tensor. Therefore, it can stabilize the algo- 
rithm of the imaginary time evolution, provided that the Trot- 
ter step r is small enough so that the bond projection oper- 
ator U(r) is nearly unitary. 27 This near unitary evolution can 
modify the wavefunction and reshape it in order to satisfy the 
canonical conditions. The simple update scheme hence pro- 
vides a quasi-canonical evolution of the iTTN state, which 
will finally converge to the ground state and become canoni- 
cal in the limit r — > 0. In practical calculations, the Trotter 
step t is gradually decreased from 10 -1 to 1CT 4 , and the total 
number of projections steps varies from 20000 to 200000. 

Now let us consider how to implement the simple update 
scheme efficiently. A simple approach is to do directly the 
singular value decomposition the evolving block tensor, which 
is a matrix of D 2 x D 2 . The computational cost for doing this 
singular value decomposition is high, since it scales as 0(D 6 ). 
This cost can in fact be reduced to 0{D A ) if we carry out this 
singular value decomposition in the following steps (again, 
projection on x bond is taken as an example): 

(1) Define the following two D 2 x Dd block matrices (Fig. 
2b) 

(At )...,:».,» = ^ylzKyv ( 4 ) 

{M b ) yXtX>m = A y7z B m xyv (5) 

by absorbing the diagonal matrices A y and y z into the tensors 
A and B, and calculate their QR decomposition 

k 

where a — a or b. Q a is a D 2 xDd column orthonormal matrix. 
R a is a Dd x Dd upper diagonal matrix. 

(2) Apply the bond projection operator U(f) to the system 
and define the gate matrix (Fig. 2c) 

G mikvJ „ 2kl = Yj < m ^ m 2Mr)\m[m' 2 )Rl.^, x e x Rl 2 . Kx (7) 

x,m\ 

(3) Take the singular value decomposition for this matrix (Fig. 
2d), 

= U mi k 1 j6'iV m2 k 1 -i, (8) 

where U and V are two Dd x Dd unitary matrices, and 6' is a 
semi-positive defined matrix. 

(4) Truncate the inner bond dimension by keeping the 
largest D matrix elements of 8', and update the local tensors 
by the formula (Fig. 2e) 

k 

Combining this efficient simple update scheme and the lo- 
cal determination of physical observables using the canonical 
form, we can keep the computational cost in a low level. In 
practice, this allows us to keep a relative large bond dimen- 
sion. 
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FIG. 3: (Color online) (a) The longitudinal and transverse magne- 
tizations m- and m x versus the transverse fields h x . The transverse 
magnetization m x increases monotonously with h x , while the lon- 
gitudinal magnetization m z decreases and vanishes at the transition 
point, (b) The second-order derivative of the ground state energy per 
site e with respect to h x , cpe/dh 2 , obtained by taking the first-order 
derivative of m x , dm jdh x . 

HI. THE TRANSVERSE ISING MODEL 

The transverse Ising model is defined by the Hamiltonian 

//,,- - ^ ./.v;.v, ^//,.v; X^'i- 

<i,j> i / 

where the spin-spin exchange constant / is set as the energy 
scale (7=1, ferromagnetic coupling). The second term repre- 
sents the transverse-field along S x direction, and the last term 
is the longitudinal-field along S : direction. 

Figure 3a shows the longitudinal and transverse magnetiza- 
tions, m, = (S z ) and m x = (S x ), as a function of the transverse 
field h x . A continuous order-disorder phase transition is found 
at h c . For h x < h c , the ground states undergo a spontaneous Z% 
symmetry breaking with a finite longitudinal magnetization 
m z , which decreases with increasing h x and vanishes at the 
critical field. By utilizing the Hellmann-Feynman theorem, 
the second-order derivative of the ground state energy can be 
calculated by d 2 e/dh 2 = -dm x /dh x . As shown in Fig. 3b, 
d 2 e/dh 2 exhibits a discontinuity at h c , indicating that h c is a 
second-order phase transition point. The critical field is found 
to be h c =2 1 .1 15, in agreement with previous calculations. 25 ' 26 
It is also close to the critical field h c = 1.06625(2) for the 
transverse Ising model on the honeycomb lattice. 28 

Figure 4 shows the bipartite entanglement entropy Se 

S E = -Tr[A 2 log 2 (A 2 )], (12) 
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FIG. 4: (Color online) (a) The entanglement entropy Se of the 
ground state for the transverse Ising model. The cusp at h c » 1.115 
corresponds to the second-order phase transition point, (b) The cor- 
relation length £ of the ground state, which also shows a cusp at the 
transition point. 



and the correlation length £ for the ground state. In Eq. (12), 
A = 8, A or y is a diagonal matrix that satisfies the canonical 
condition. S E shown in Fig. 4a is obtained by taking the 
average over the three bonds. 

The correlation length £ is evaluated from the ratio of the 
the largest (do) and the second largest eigenvalue (a\) of the 
transfer-matrix for the iTTN state, 



f=l/ln2» 



(13) 



For the q — 3 Bethe lattice, there are six kinds of transfer- 
matrices, depending on the site and the bond direction. For 
example the transfer matrix along the yz-direction is defined 
by (for A sublattice site) 



fit 



The other five transfer matrices including r* 



pa,b 



T"' rx , y can be similarly defined. The results of the correlation 
length shown in Fig. 4b are evaluated from the product of the 
six transfer matrices along a specific path. 

A distinctive feature revealed by Fig. 4 is that the correla- 
tion length as well as the entanglement S e, does not diverge 
at the critical point. £ is found to be upper bounded by 1 / In 2, 
in agreement with the published results 25,26 . This peculiar be- 
havior is not observed in the ordinary continuous phase transi- 
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FIG. 5: (Color online) The field dependence of the longitudinal sus- 
ceptibility x z ■ It shows a divergent peak at the transition point h c ^ 
1.116. The susceptibility is calculated by x z = [> n z(hz) ~ m(0)]/Sh 
With h, = 10- 4 . 



tion systems, where the correlation length is always divergent 
at the critical point. 

We will now show that this non-critical behavior of the cor- 
relation length at the critical point is due to the peculiar ge- 
ometry of the Bethe lattice. For the Bethe lattice, the num- 
ber of sites on the boundary of a finite connected region is 
roughly equal to the number of internal sites within that re- 
gion. It means that the lattice sites are highly non-uniformly 
distributed as a function of lattice distance away from a given 
center. This is a feature of the Bethe lattice that differs from a 
regular lattice. 

In order to understand why the correlation length is finite at 
the critical point, let us take a scaling transformation to con- 
vert the Bethe lattice to a "regular" 2D lattice whose lattice 
sites are uniformly distributed in space. To do this, we first 
choose an arbitrary site, to be viewed as "center" of the lat- 
tice, and define the distance R for a given layer r to the center 
as 



R oc 



N(r) 



(q-lY 



(15) 



where N(r) oc (q - 1)'' is the number of sites enclosed by the 
r layer. In this rescaled lattice, r ~ 2 InR + const, and the ex- 
ponentially decaying correlation function C(r) ~ exp(-r/£) in 
the original Bethe lattice corresponds to a power-law decaying 
function of R 



C(R) ~ R 



(16) 



The algebraic decay of this correlation function suggests the 
spins on the Bethe lattice are actually long range correlated 
in the rescaled framework, even away from the critical point. 
This is the reason why the system can undergo a phase tran- 
sition without exhibiting a divergent correlation length at the 
critical point in the original Bethe lattice. 
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To understand why the correlation length is upper bounded, 
let us consider the longitudinal magnetic susceptibility Xz of 
the ground state (Fig. 5) 



GS 



<^>GS<Sr> 



GS- 



(17) 



where S i is the spin at a reference center, i runs over all the 
sites on the lattice. (O)cs is the expectation value of opera- 
tor O in the ground state. Exploiting the Cj, rotational sym- 
metry of the Bethe lattice, we define C z (r) = (SiS z )gs - 
(S^)gs(S])gs, with r the layer number where site i resides. 
Thus Xz can t> e rewritten as 



Xz = 2 n(r)C z (r), 



(18) 



where n(r) is the number of spins on layer r. On a regu- 
lar lattice, n(f) oc r v ~', where v is the spatial dimension of 
the lattice, the susceptibility is always finite if the spin-spin 
correlation function C z (r) decays exponentially. However, 
in the Bethe lattice, n(r) oc (q - l) r_1 . Now if we assume 
C z (r) oc exp(-r/£), then 



^«J](9-1)V^ = 2' 



r[ln(?-l)-l/f] 



(19) 



which diverges if £ approaches the threshold value 1 / ln(^- 1). 
This shows that the susceptibility can diverge even if C z (r) 
decays exponentially with r. The critical point occurs when 
£ = 1/ ln(^- 1), and the correlation length £ is therefore upper- 
bounded by 1/ \n(q - 1) on the Bethe lattice. 29 



IV. ANISOTROPIC HEISENBERG MODEL 

The anistropic Heisenberg model, i.e. the XXZ model, is 
defined by the Hamiltonian 



(20) 



<i,j> 



where 6 is the anisotropy parameter. 

The above model has been intensively studied on the hon- 
eycomb and square lattices by different numerical meth- 
ods, which include exact diagonalization, 30 quantum Monte 
Carlo, 30,31 coupled cluster methods, 32 ' 33 and tensor network 
algorithms. 12,34 It is found that the system possesses magnetic 
long-range orders for all values of 6. The antiferromagnetic 
ordering vector points within the easy xy-plane for 6 < 1 or 
along the z-axis for 5 > 1 . There is a first-order phase transi- 
tion at 5 — 1, the Heisenberg point. 35 

This model was also studied on the Bethe lattice (more pre- 
cisely, on the Cayley tree lattice) by DMRG, 21 23 ' 36 . It was 
found that there exists a long-range magnetic order at the 
isotropic point 5 — 1. It was also suggested that a quantum 
phase transition occurs at this point. However, the properties 
of this transition and the phases on the two sides of the critical 
point have not been clarified. 
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FIG. 6: (Color online) (a) The ground state energy and (b) the bipar- 
tite entanglement entropy as a function of d for the XXZ model. The 
appearance of the cusp in the energy.as well as the discontinuity in 
the first-order energy derivative [inset in (a)], suggests that this is a 
first-order phase transition point. 
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FIG. 7: (Color online) The staggered magnetization along the z-axis 
(mi) and x-axis (ml.). The spin orientations flip suddenly at the transi- 
tion point 5=1. For 5 < 1 , the spins are ordered within the xy-plain, 
while for 6 > 1, the ordering is along the z-axis. 



Figure 6 shows the ^-dependence of the ground state energy 
per bond and the entanglement entropy for the XXZ model. A 
clear first order quantum phase transition is observed at 6 — 1 . 
The energy per bond shows a change of slope at 6 = 1 (the 
first-order energy derivative is shown in the inset), which sug- 
gests that there is an energy level crossing. The entanglement 
entropy varies continuously across the transition point, but ex- 
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FIG. 8: (Color online) The one-ring cluster with 12 sites. The diag- 
onal matrices A are defined on the dangling bonds. The dashed lines 
represent removed bonds (the physical couplings on the dashed line 
still exist). 



FIG. 9: (Color online) The 4-ring cluster with 26 sites. The inequiv- 
alent lattice sites are numbered from 1 to 16. Ay labels the diagonal 
matrix on the bond linking sites i and j. 



hibits a cusp. 

Figure 7 shows the staggered magnetization m\ and mi 
around the critical point. The ground state is found to possess 
in-plane antiferromagnetic order with a finite m\ for 6 < 1, 
and z-axis antiferromagnetic order with a finite mi for 6 > 1 . 
At the transition point, the two order parameters change sud- 
denly, displaying a spin flip transition. This result verifies the 
conjecture made by Otsuka. 21 

At the transition point, we find that the ground state en- 
ergy per bond has the value e/, = -0.359817(3) and the spon- 
taneous magnetization has the value m s = 0.34736(1) for 
D = 40. The errors in the parentheses are estimated by com- 
paring the results for different bond dimension D and differ- 
ent Trotter slices r. Our results agree well with the DMRG 
data published in Ref. [23], where the local magnetization is 
found to be m = 0.347 on the central lattice site and the bond 
energy between the central spin and a spin on the first layer 
is e — -0.359. This satisfactory agreement suggests that by 
calculating the Bethe lattice, we can reproduce the results of 
local properties in the very center of a large Cayley tree. 



V. CLUSTER UPDATE SCHEME 

In the previous sections, the simple update has been ap- 
plied to study the quantum spin models on the Bethe lattice, 
leading to very accurate results. What is more, in terms of 
the Bethe approximation, these results can also be regarded 
as approximations for the corresponding 2D lattice models. 
Actually, the simple update scheme has already been used to 
study regular 2D lattices, such as the honeycomb or square lat- 
tices. Combined with the TRG/SRG techniques, it can achieve 
rather accurate results. 7 - 910 Nevertheless, in this section, we 
would provide a different way of using the simple update to 
calculate 2D lattices. Inspired by the generalization of the 




— e — Simple update, D=40 
QMC 
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FIG. 10: (Color online) The ground state energy per site for the 
Heisenberg model on the honeycomb lattice. The results are eval- 
uated in the central area of the clusters. The dashed line is the re- 
cent quantum Monte Carlo result. 37 The number of hexagonal rings 
is used to label the cluster size. 



Bethe approximation to larger clusters in classical statistical 
mechanics, 17 we apply simple update to various tree tensor 
clusters. In this way, the advantage (efficiency) in treating 
a tree tensor network, namely the fact that it can be readily 
canonicalized, is utilized to improve the calculation accuracy 
on a regular 2D lattice. 

To start, as a first-order approximation, let us compare the 
results on the q — 3 Bethe lattice (which has no loops) with 
those on the 2D honeycomb lattice (whose coordination num- 
ber is also q — 3, and it does have loops). Our result for the 
ground state energy of the Heisenberg model on the q — 3 
Bethe lattice is e\, = -0.359817(3), while the correspond- 
ing energy on the honeycomb lattice obtained by the recent 
quantum Monte Carlo calculation is cqmc = -0.36303(14). 37 
The relative difference between these two energies is less 



8 



than 0.9%. However, the spontaneous magnetization for the 
ground state of the Heisenberg model on the Bethe lattice, 
namely m s = 0.34736(1), is much larger than the correspond- 
ing value on the honeycomb lattice, which is about 0.27 as ob- 
tained by the quantum Monte Carlo. 37 Notice, some other re- 
sults for the magnetization on the honeycomb or square lattice 
obtained with the tensor network algorithms are also found to 
be higher than the Monte Carlo ones. 7 ' 34 

As a next step, our approximate treatment of the 2D hon- 
eycomb lattice can be improved by using tensor networks that 
include rings. In Fig. 8, the cluster with one hexagonal ring 
is shown, some geometric bonds are removed (dashed lines in 
Fig. 8) to form a tree tensor cluster. Note although the tensor 
network does not have geometric bonds on the dashed lines, in 
the Hamiltonian the couplings along these bonds nevertheless 
exist. Therefore, the projections by imaginary time evolution 
should be executed also on the dashed lines. This cannot be 
done directly as on usual bonds, but can be accomplished as 
follows with the help of the swap gates. The swap gates are 
used to exchange the physical indices of two tensors, which 
proceeds similarly as the projection scheme illustrated in Fig. 
2, with the minor revision that the imaginary tim evolving op- 
erator U (t) is now replaced with a swap operator U s , that con- 
ducts U s \nii,mj) — \mj,m.i). 

In Fig. 8, take the dashed bond between site A and F as 
an example, swap gates moves the physical index on site A 
in the order A — > B — > C, and the physical index on site F as 
F — > E — > D. After that, the two spins are linked by the solid 
bond between C and D, then we can take the projection and 
update processes as on an usual bond. After that, we have 
to move the two spin indices back to their original positions 
by reversed swap operations, which accomplishes the special 
projection step on a dashed bond. Through iterative and self- 
consistent projection processes on the solid and the dashed 
bonds of the tree cluster, an approximation for 2D lattices can 
be obtained. Compared with the simple Bethe lattice, this tree 
tensor cluster approach can provide better approximation for 
2D. On the other hand, it can also be regarded as an ideal 
method for evaluating the "super Bethe lattice", of which each 
"single site" is now placed with a hexagonal ring, instead of a 
single site, and the coordinate number q — 6. 

Beyond the one-ring cluster, more rings can be included 
to further enlarge the cluster. As an example, Fig. 9 shows 
a cluster with 4 hexagonal rings. The accuracy of energy 
versus different cluster size (labelled by the number of rings 
included) are shown in Fig. 10, which verifies that the ac- 
curacy could be improved consistently with enhancing the 



cluster size. To obtain better approximation for true 2D lat- 
tices, the local observables are detected in the center area 
of the cluster. In practice, for the 4-ring tree tensor clus- 
ter in Fig. 9, the results are obtained by averaging over 
sites 3 and 4. We find an energy per site of e -0.54441 
(bond energy e/, =* -0.36294) and a local magnetization of 
m - [e(h s ) - e(h s = Q)]/h s - 0.3147 (with a staggered mag- 
netic field h s = 0.01). Hence, the inclusion of rings clearly 
improves the agreement with QMC data. For the transverse 
Ising model, through the 4-ring cluster calculations, the phase 
transition point is estimated as h c = 1.1, which is also more 
accurate than the simple Bethe approximation. More numeri- 
cal results with larger clusters and further details of this cluster 
Bethe approximation will be published separately. 



VI. CONCLUSION 

In summary, the simple update scheme is employed to study 
two spin models on the Bethe lattice, i.e., the transverse Ising 
and the Heisenberg XXZ model. For the Ising model, it is 
found that the correlation length, as well as the entanglement 
entropy, does not diverge at the second-order transition point. 
Through a scale transformation, we have given an intuitive ex- 
planation of this peculiar "critical" phenomenon. Moreover, 
by studying the magnetic susceptibility, we show that the cor- 
relation length is upper bounded. For the Heisenberg XXZ 
model, the existence of a first-order phase transition at the 
isotropic point is clearly verified, and the two different mag- 
netic ordered phases are identified as the easy-plain and easy- 
axis phases, respectively. Furthermore, in terms of the Bethe 
approximation, we obtain accurate and scalable approxima- 
tions for the 2D lattice models by applying the simple update 
to tree tensor clusters. 
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